w t y l j r r^ * 


Discretization and Preconditioning Algorithms for 
the Euler and Navier-Stokes Equations on 
Unstructured Meshes 


Tim Barth 
NASA Ames R.C. 
Moffett Field, California USA 


Contents 


1 Symmetrization of Systems of Conservation Laws 3 

1.1 A Brief Review of Entropy Symmetrization Theory 4 

1.2 Symmetrization and Eigenvector Scaling 5 

1.3 Example: Compressible Navier-Stokes Equations 6 

1.3.1 Entropy Scaled Eigenvectors for the 3-D Euler Equations 8 

1.4 Example: Quasi- Conservative MHD Equations 9 

1.4.1 Powell’s Quasi-Conservative Formulation of the MHD Equations . . 10 

1.4.2 Symmetrization of Powell’s Quasi-Conservative MHD Equations . . 11 

1.4.3 Eigensystem of Powell’s Quasi-Conservative MHD Equations .... 12 

1.4.4 Entropy Scaled Eigensystem of Powell’s Quasi-Conservative MHD 

Equations 12 

2 Maximum Principles for Numerical Discretizations on Triangulated Do- 
mains 

2.1 Discrete Maximum Principles for Elliptic Equations 16 

2.1.1 Laplace’s Equation on Structured Meshes 16 

2.1.2 Monotone Matrices 17 

2.1.3 Laplace’s Equation on Unstructured Meshes 17 

2.2 Discrete Total Variation and Maximum Principles for Hyperbolic Equations 21 

2.2.1 Maximum Principles and Monotonicity Preserving Schemes on Mul- 
tidimensional Structured Meshes 25 

2.2.2 Maximum Principles and Monotonicity Preserving Schemes on Un- 
structured Meshes 26 

3 Upwind Finite Volume Schemes 34 

3.1 Reconstruction Schemes for Upwind Finite Volume Schemes 34 

3.1.1 Green-Gauss Reconstruction 35 

3.1.2 Linear Least-Squares Reconstruction 36 

3.1.3 Monotonicity Enforcement 38 

3.2 Numerical Solution of the Euler and Navier-Stokes Equations Using Upwind 

Finite Volume Approximation 42 

3.2.1 Extension of Scalar Advection Schemes to Systems of Equations . . 42 

3.2.2 Example: Supersonic Oblique Shock Reflections 43 

3.2.3 Example: Transonic Airfoil Flow 44 

3.2.4 Example: Navier-Stokes Flow with Turbulence 44 


i 



4 Simplified GLS in Symmetric Form 48 

4.1 Congruence Approximation 48 

4.2 Simplified Galerkin Least-Squares in Symmetric Form 49 

4.2.1 A Simplified Least-Squares Operator in Symmetric Form 49 

4.2.2 Example: MHD Flow for Perturbed Prandtl-Meyer Fan 50 


2 



Chapter 1 

Symmetrization of Systems of 
Conservation Laws 


This lecture briefly reviews several related topics associated with the symmetrization of 
systems of conservation laws and quasi-conservation laws: 

1. Basic Entropy Symmetrization Theory 

2. Symmetrization and eigenvector scaling 

3. Symmetrization of the compressible Navier-Stokes equations 

4. Symmetrization of the quasi-conservative form of the MHD equations 

There are many motivations, some theoretical, some practical, for recasting conservation 
law equations into symmetric form. Three motivations are listed below. The first motivation 
is widely recognized while the remaining two are less often appreciated. 

1. Energy Considerations. Consider the compressible Navier-Stokes equations in quasi- 
linear form with u the vector of conserved variables, f* the flux vectors, and M the 
viscosity matrix 

u ,t + f u u ,s* = (MfjU X j (^*1) 

In this form, the inviscid coefficient matrices F u are not symmetric and the viscosity 
matrix M is neither symmetric nor positive semi-definite. This makes energy analysis 
almost impossible. When recast in symmetric form, the inviscid coefficient matrices 
are symmetric and the viscosity matrix is symmetric positive semi-definite. The energy 
analysis associated with Friedrichs systems of this type is well-known. 

2. Dimensional consistency. As a representative example, consider the time derivative 
term from (1.1). The weak variational statement associated with this equation requires 
the integration of terms such as — / w^u dxdt. When w and u reside in the same 
space of functions, the inner product quantity w T u is dimensionally inconsistent. 
Consequently, errors made in a computation would depend fundamentally on how 
the equations have been non-dimensionalized. When recast in symmetric form, the 
inner product w r u is dimensionally consistent with units of entropy density per unit 
volume. 
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3. Eigenvector scaling. Apart from degenerate scalings, any scaling of eigenvectors satis- 
fies the eigenvalue problem. Unfortunately, numerical discretization techniques some- 
times place additional demands on the form of right eigenvectors. As noted by Balsara 
[Bal94] in his study of high order Godunov methods, several of the schemes he stud- 
ied that interpolate “characteristic” data (see for example Harten et al. [HOEC87]) 
showed accuracy degradation that depended on the specific scaling of the eigenvectors. 
In the characteristic interpolation approach, the solution data is projected onto the lo- 
cal right eigenvectors of the flux Jacobian, interpolated between cells, and finally trans- 
formed back. The interpolant thus depends on the eigenvector form. Symmetrization 
provides some additional insight into the scaling of eigenvectors. Let R(n) denote the 
matrix of right eigenvectors associated with the generalized flux Jacobian in the direc- 
tion n. Using results from entropy symmetrization theory, Sec. 1.2 describes a scaling 
of right eigenvectors such that the product R(n)R T (n) is independent of the vector n. 
This result is used in Sec. 1.4 which discusses the ideal magnetohydrodynamic (MHD) 
equations. The right eigenvectors associated with these equations exhibit notoriously 
poor scaling properties, especially near a triple umbilic point where fast, slow and 
Alfen wave speeds coincide, [BW88]. The entropy symmetrization scaling provides 
a systematic approach to scaling eigenvectors which is unique in the sense described 
above. 

1.1 A Brief Review of Entropy Symmetrization Theory 

Consider a system of m coupled first order differential equations in d space coordinates and 
time which represents a conservation law process. Let u(x,t) : H d x R + i-> R m denote the 
dependent solution variables and f(u) : R m »-* R mXd the flux vector 

u ,£ + = 0 ( 1 - 2 ) 

with implied summation on the index i. Additionally, this system is assumed to possess the 
following properties: 

1. Hyperbolicity. The linear combination 

f,u(n) = f* u 

has m real eigenvalues and a complete set of eigenvalues for all n € R d . 

2. Entropy Inequality. Existence of a convex entropy pair U(vl),F(u) : R m > R such 
that in addition to (1.2) the following inequality holds 

U,t + F: Xi < 0. (1.3) 

In the standard symmetrization theory [God61, Moc80, Har83b], one seeks a change of 
variables u(v) : R m R m to Eqn. (1.2) so that when transformed 

u, V v, t + r v v, Xj = 0 (1.4) 

the matrix u iV is symmetric positive-definite and the matrices f' v are symmetric. This 
would be a classical Friedrichs system. Clearly, if functions £f(v),.P(v) : R" 1 i-> R l can be 
found such that 

u = u, y , r - .7> 
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then the matrices 


U ?v — ^,V,V) fv ~~ ^v,v 

are symmetric. To insure positive-definiteness of u jV so that mappings are one-to-one, 
convexity of U{\) is imposed. Since v is not yet known, little progress has been made but 
introducing the Legendre transform 

U( u) = u r v — U(\) 

followed by differentiation 

^ u = V T +U T V, u -W v V u =V T 

yields an explicit expression for the entropy variables v in terms of the entropy function 
U(u). Symmetrization and generalized entropy functions are intimately linked via the 
following two theorems: 

Theorem 1.1.1 Godunov [God61] If a hyperbolic system is symmetrized via change of 
variables, then there exists a generalized entropy pair for the system. 

Theorem 1.1*2 Mock [Moc80] If a hyperbolic system is equipped with a generalized en- 
tropy pair f/, then the system is symmetrized under the change of variables v T = U }U . 

For many physical systems, entropy inequalities of the form (1.3) can be derived by ap- 
pealing directly to the conservation law system and the second law of thermodynamics. 
Using this strategy, specific entropy functions for the Navier-Stokes and MHD equations 
are considered in Secs. 1.3, 1.4 respectively. 

1.2 Symmetrization and Eigenvector Scaling 

In this section, an important property of right (or left) symmetrizable systems is given. 
Simplifying upon the previous notation, let A 0 = u jV ,vl l = f* v and rewrite (1.4) 

^,,+^,=0. (1.5) 

SPD Symm 

The following theorem states an important property of the symmetric matrix products A 1 A 0 
symmetrized via the symmetric positive definite matrix .4° . 

Theorem 1.2.1 (Eigenvector Scaling) Let A £ R” * " be an arbitrary diagonalizable 
matrix and S the set of all right symmetrizers: 

S = {B € R" xn | B SPD, AB symmetric}. 

Further, let R £ R/ £Xn denote the right eigenvector matrix which diagonalizes A 

A = RAR 1 

with r distinct eigenvalues, A = Diag(Ai , Wm-ixm-ii ■ • • , KIm r xm r )- Then for each 

B £ S there exists a symmetric block diagonal matrix T = Diag(T m , xm , , T irt2 xm ., , . . . , T 1UrXmr ) 
that block scales columns of R. R = R'I\ such that 

B = RR r , A = RAR 1 


b 


which imply 


AB = RAR 1 '. 


Proof: The symmetry of B and A B implies that 

AB-BA t = RAR-'B - BR~ t AR t = 0 
or equivalently for Y € R nxn 


AT - FA = 0, Y = R~ l BR T . (1.6) 

Partition Y into r x r blocks, Y miX m^ with block dimensions corresponding to eigenvalue 
multiplicities. Equation (1.6) then reduces to the following set of decoupled systems: 

Imi x mi x nij ^jYimxmj Imj xmj — 0? — T. (1*7) 

or simply 

(Xi — Xj)Y mi xmj = 0, (1.8) 

This implies that Y is of block diagonal form since Y 1JliXm . = 0,i ^ j. Prom the definition 
Y — R~ l BR~ T , Y is congruent to hence symmetric positive definite (SPD). Given the 
block diagonal structure of 7, the square root factorization exists globally as well as for each 
diagonal block Y m%xmi = This y* e lds the stated theorem with T = T 1 / 2 . 

■ 

This theorem is a variant of the well-known theory developed for the commuting matrix 
equation 

AX-XA = 0, A,X€R nxn , 

see for example Gantmacher [Gan59], Note that the general theory addresses the more 
general situation for which the matrix A can only be reduced to Jordan canonical form. 


Remark 1.2.1 Prom the scaling theorem 1.2.1, the right eigenvectors associated with each 
A 1 can be scaled so that A 1 A 0 = R l Ai(R l ) T which yields a revealing form of the symmetric 
quasilinear form 


A\ t + R i A i (R i ) 1 v, Xi = 0 
with A 0 = R 1 {R l ) T = • • • = R d {R d ) T . 


1.3 Example: Compressible Navier-Stokes Equations 


Consider the ideal compressible Navier-Stokes equations, x € R d , 


d_ 

dt 



+ V 


pV 

p\\ -hip 
V{E + p) 


V- 


rV-q, 


(1.9) 


where V € R d is the velocity vector, p and p the density and pressure of the fluid, E the 
specific total energy defined as 


E = 


P 

7 - 1 


+ J/’V 2 


( 1 . 10 ) 


6 



and t is the viscous stress tensor: 


T = A 





( 1 . 11 ) 


In addition, an ideal gas is assumed p = p RT as well as Fourier heat conduction q = — k VT. 
In these equations A and p are diffusion coefficients, k the coefficient of thermal conductivity, 
7 the ratio of specific heats, and R the ideal gas law constant. 

In 1983, Harten [Har83b] proposed the generalized convex entropy function 


U(u) = - pg{s ), 


, q 
9 > 0 , — < 7 

9 


l 


for the compressible Euler equations. In this equation, s is the thermodynamic entropy of 
the fluid. This choice was motivated from the well-known entropy transport inequality for 
the inviscid Euler equations: 

s y t + V • Vs > 0 


which generalizes to 


g(s),t + V • Vg(s) >0, g' > 0. 


or after combining with the continuity equations 


(p9( 3 )),t + V • pg{s)V > 0. 

Thus, it becomes clear that U = —pg(s),F l = -pg{s)Vi is an exceptable entropy pair. 
Hughes, Franca, and Mallet [HFM86] removed the arbitrariness of g(s) by showing that 
symmetrization of the Navier-Stokes equations with heat conduction places the additional 
restriction that g(s) be at most affine in s, i.e. g(s) = cq + c\s. A convenient choice is given 
by U(s) = which yields the following entropy valuables: 

/ _l_ 2±i _ K \ 

7-1 ' 7 — 1 P \ 

v = U,u= f 

v -% > 

The change of variable matrix u jV takes a particularly simple form: 

(9 pV T E \ 

u iV = pV pYV + pi pH\ 

\E pHV T pH 

with a the sound speed, a 2 = 7 p/p, and H the specific total enthalpy, H — a 1 j { 7 — l)+V 2 /2. 

Consider the application of the Scaling Theorem 1.2.1 to the inviscid Euler terms ap- 
pearing in the Navier-Stokes equations 

A°v it + A i A°v iXi = 0 ( 1 . 12 ) 

with A 1 = f* u and A 0 = u v . It is sufficient to consider symmetrization of the arbitrary 
linear combinations of the form 


d(n) = niA\ ||n|| = 1 


(1.13) 
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and scaled right eigenvectors R( n) such that 

A(n) = R(n)A(n)R- l (n), A 0 = R(n)R r (n). 

From this result, the symmetric coefficient matrices are given by 

A(n)A° = R(n)A(n)R T (n). 

Prom a derivational point-of-view, it is advantageous to first compute the eigensystem as- 
sociated with the system in primitive variables, w = (p, V,p) r , 

w,t + Uw A 1 u w w tXl = 0 (1.14) 

and then to compute the scaling of the right eigenvectors, r(n), of the primitive variable 
system. Once the scaled right eigenvectors of the primitive variable system have been 
computed, scaled right eigenvectors of the conservative system are easily recovered from 


with 


R( n) = 

u, w r\ 

( 1 

0 T 

u, w = V 

pi 

\Jv* 

pV T 



(1.15) 


1.3.1 Entropy Scaled Eigenvectors for the 3-D Euler Equations 

Following the procedure described in Theorem 1.2.1, after some tedious algebraic manip- 
ulation, the following scaled right eigenvectors of the primitive variable system have been 
obtained: 


Entropy and Shear Waves: Ai 2,3 = V * n 


n,2,3 = 



y/y-lpn 

a[C] 

0 T 



where [C(n)] = n^e and is the usual alternation tensor. 
Acoustic Waves: A± = V • n ± a 


(1.16) 


r± 



(1.17) 


In Fig. 1.1, the FAiclidean norm, ||/?(n)7? T (n)||£, is graphed for n = (cos 0, sin0 ) r , 8 £ 
[0,27t] with constant (p, V,p) using entropy scaled eigenvectors and the “naturally” scaled 
eigenvectors given in Struijs [Str94]. The entropy scaled eigenvectors produce a constant 
result since R(n)R r (n) — u jV which is independent of n. Although the naturally scaled 
eigenvectors seem to differ only in minor ways from the entropy scaled counterparts, the 
MHD example given in the next section provides a more convincing argument for proper 
eigenvector scaling. 
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Figure 1.1: Euclidean norm dependency of R(n)R(n) T on n, {pjrho 0 0 = 1 , Vj / = 

•8,^2/000 = -.6, p/poo = 1.4) 


1.4 Example: Quasi-Conservative MHD Equations 

As a second example, consider the ideal (nonrelativistic) MHD equations, x £ R d , 


l P \ 

d_ pV 
dt E I 

V B / 


pV \ 

pVV + I (p + |B 2 ) - BB 

V (E + p + ±B 2 ) -B (V B) 

VB - BY / 


(1.18) 


In addition to the variables defined for the Navier-Stokes equations, B E R d is the magnetic 
field, and E is the specific total energy redefined as 

E=-^r+p\V> + \ B 2 . (1.19) 

7—1 Z L 

At this point, the equations have been written in conservative (divergence) form. Unlike 
the Navier-Stokes equations, the MHD equations have the additional time independent 
constraint V • B = 0. A number of numerical techniques exist for solving problems of this 
type: 


L Mixed finite element formulations and Lagrange multiplier formulations 


2. Finite element penalty formulations. 

3. Staggered mesh formulations. 

4. Projection formulations. 

In the next section, a technique developed by [Pow94] is used to reformulate the MHD 
equations in quasi-conservative form which obviates the V • B = 0 constraint. 

In has been known for some time that negated specific entropy, — ps, is a convex gener- 
alized entropy function for the ideal MHD equations, see Rugger i [RS81]. Even so, starting 
from the MHD equations (1.18) in the standard quasi- linear form: 


U t + /l'u z, = 0, 


yl l = f* 

*,11 


( 1 . 20 ) 
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followed by the change of variables to 


A°v tt + ^A\ Xi = 0, A° = u v (1.21) 

with v r = U, u and f/(u) = — does not symmetrize the system. For example in 2-D, the 


following result 

is obtained: 

o 

0 

0 

0 

0 

0 

\ 



0 

0 

0 

pB? 

p 

pB 1 

p 

0 




0 

0 

0 

pB\ Z?2 

pB 2 

0 


- 

(A l A°f = 

0 

Pil 

pBi B 2 

p 

0 

P 

pV 2 B 2 

pV 2 B\ 




p 

p 

p 

P 




0 

pfi 1 

pB 2 

pV 2 B 2 

0 

pV 2 




p 

p 

p 

P 




Vo 

0 

0 

pViBi 

p 

p 

0 

) 


This failure to symmetrize the MHD system can be explained by the simple observation 
that the constraint V-B — 0 has not been used. A straightforward derivation of the entropy 
transport equation for the ideal MHD equations reveals this. For ideal MHD, entropy is 
given by s = log(pp 7 ) so that 

7 1 

ds = dp + -dp. 

P P 

Inserting terms from (1.18), the following transport equation for smooth flow results: 

s jt + v • Vs + (7 - l)(v • B)(V • B) = 0. (1.23) 

Clearly, the V * B = 0 condition is fundamental in the derivation of the generalized entropy 
function. The next section considers Powell’s method which implicitly satisfies V * B — 0 
and permits symmetrization using U (u) = — 

1.4.1 Powell’s Quasi-Conservative Formulation of the MHD Equations 

During an investigation of approximate Riemann solvers for the ideal MHD equations, 
Powell [Pow94] observed that the coefficient matrices, A\ appearing in the quasi-linear form 
are individually rank deficient, i.e. each contains a zero row corresponding to components of 
the B field. This degeneracy was removed by replacing the zero eigenvalue and eigenvector 
with a “divergence wave” eigenvector with velocity eigenvalue similar to the entropy wave. 
Powell then noted that adding the divergence wave eigenvector was equivalent to writing 
the MHD equations (1.18) in chain rule form: 


dpV 

~dt~ 


! + v.(„ v )=o 


o 


+ V • (pW) + V (p + ^B 2 ) - B VB - B(V-B) = 
^+V- v(tf + p + l B 2 ) + (e + p + ^B 2 ) V V - B V(V • B) - (V ■ B)(V • B) = 0 


m 

dt 


4-VVB+BVV-B VV- V( V ■ B) = 0 


and weakly enforcing V * B — 0 by removing terms proportional to V ■ B as underlined in 
Eqn. 1.24. Powell also rioted that this modification changes the nature of the equations in 


(1.24) 
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a fundamental way. This is revealed by taking the divergence of the B field equations for 
the original system: 

V ' (^ + V (VB_BV) ) = f (V ‘ B)= ° (L25) 

as well as the modified system (1.24) with underlined terms removed: 

vf? + VVB + BVV-BVv'| = ^(V-B)+V-(VV-B) =0. (1.26) 

\ at /at 

The first form (1.25) states that if V • B = 0 is initially zero then it should remain zero. 
The second form (1.26) states that (V • B)/p is a passive scalar for the system. Any local 
V-B is simply advected away. Powell asserts that this is a more numerically more stable 
process. Finally, note that removal of the underlined terms in (1.24) can be viewed as 
adding source-like terms proportional to V • B to the original divergence form (1.18), i.e. 



+ V 


( pV 

pVV + 1 (p+ iB 2 ) - BB 

V (e + p+ iB 2 ) — B (V-B) 
V VB - BV 



0 \ 

B 

V-B 


V / 


(V • B). (1.27) 


This suggests the quasi-conservative nature of the equations. In the next section, the 

7-1 


coefficient matrix form is retained in order to show that the entropy function U(n) = — 


does symmetrize Powell’s modified MHD equations. 

1.4.2 Symmetrization of Powell’s Quasi-Conservative MHD Equations 

It is again most convenient to write the equations in the primitive variable form w = 
(p,V,p,B) T : 

(1.28) 


w, t + u w L A£u W w x< = 0 

where A p denotes a matrix for the modified system. Consider arbitrary combinations, 
A p ( n) = rti A l p and write Powell’s modified coefficient matrices in the following compact 
form: 

\ 


u ,wMp( n ) u ,w 


with 


/V -n 

pn T 

0 

o r 

0 


(V • n)7 

in * 

p f 

;(nB r - B n) 

0 


7pn T 

V n 

0 

0 

Bn r - B 

n 0 

(V • n)7 


1 

0 T 

0 0 7 \ 



V 

P I 

O 

o 

o 


w ~ 

1 V 2 

2 v 

1 pV T 

-i-T B T 

7-1 



0 

00 T 

0 I ) 


(1.29) 


A straightforward calculation reveals that U (u) = 
entropy variables are given by 


(1.30) 


4 does symmetrize this system. The 


V = U n = 


/ s_ + 2±i _ K + Bi\ 

/ 7 -1 ' 7 1 p ' 2 p ' 

pV 

P 

-E 
V 
pQ 

P 
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and the Riemannian matrix u v 

can be written simply as 



P 

pY T 

E-\ B 2 

0 T \ 


pV 

pYY + pi 

pHY 

00 T 

U,v = 

B 

pHY T 

pH 2 - ^ ^ 

£ B r 

P 


\ o 

O 

o 


v> 


with a the sound speed, a 2 = 7 p/p, and H the specific total enthalpy, H = a 2 j( 7— l)+V 2 /2. 


1,4,3 Eigensystem of Powell’s Quasi-Conservative MHD Equations 

In order to give the eigensystem for the MHD equations, it is useful to define b = B j yfp 
and the fast and slow speeds: 

C /,S = \ (“ 2 + b2 ) ± \\j (“ 2 + b 2 ) 2 - 4a 2 (b • n) 2 (1.31) 

The eigenvalues and eigenvectors are then written compactly as: 

Entropy and Divergence Waves: Ai 7 2 = V * n 


0 

ri= 0 ’ 

\ 0 / 

Alfv^n Waves: \± a = V • n ± (b • n) 

( 0 \ 

±(n x B) 

r±a= 0 

Vvp(nxB)/ 



(1.32) 


(1.33) 


Magneto-acoustic Waves: A ±/,± s = V • n ± c f yS 


r ±f,±s 


( 


\ 


±c f,s 


<7,s 


p 

Cf , n-(h-n) b 
*» -(b-n)» 


2 

P a 

(Bn-nB) n 

<7 ,,-( b -") 2 


\ 

/ 


(1.34) 


In this form, the magneto-acoustic eigenvectors exhibit several forms of degeneracy as care- 
fully described in Roe and Balsara [RB96]. In the next section, the entropy scaled eigenvec- 
tors are given. These eigenvectors are similar (but not identical) to the eigenvectors given 
in Roe and Balsara. The behavior of the new eigenvectors is significantly improved. 


1.4.4 Entropy Scaled Eigensystem of Powell’s Quasi-Conservative MHD 
Equations 


After consider algebraic manipulation, entropy scaled eigenvectors corresponding to the 
Powell’s quasi-conservative MHD equations have been obtained. Using the notation of Roe 
and Balsara define 


«/ = 


a 2 - C f~ n 

r.2 _ ,:i ^ - r 2 


7 


7 


(1.35) 
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and n J , a unit vector orthogonal to n lying in the plane spanned by n and b. 

Entropy and Divergence Waves: Ai 2 — V • n 


r i = 


7-1 


7 


/v/P\ 

0 

0 

V o / 


r 2 = W- 


l 0 

0 

0 

\an> 


(1.36) 


Alfv4n Waves: A± a = V • n ± b • n 



/ 0 

\ 

/T 

T*f(n L 

x n) 

r±a= \ 2 

0 


v \/! (n± 

x “) ) 

Waves: A±y = 

V * n ± cj 



r± f - 


( , a /v/P 

. a/ a 2 n+Qj o((b-n i )n-(b-n)n 1 ) 

vW 
a f ,/pa 2 

\ ajon 1 


(1.37) 


(1.38) 


Slow Magneto-acoustic Waves: X± s = V • n ± c, 

/ <*sy/p 

nr i -J-semfh * n) 

r± s = 


\ 


. a 3 a(b-n) n+a/ cj n 1 - 

±sgn(b-n) ^ c/ 7 


a. 


V 


VP a 


“ a / a 


n" 


/ 


(1.39) 


Next, the experiment performed in Sec. 1.3.1 is repeated for the MHD equations. In Fig. 



Figure 1.2: Euclidean norm dependency of R(n)R(n) 7 on n for the MHD equations, 
(p/rhooo = 1.3, Vi/aoo — .8,V 2 / aoo = — -^iP/Poo = 2.0, R x = ,2,I3y = 1.2) 
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1.2, the Euclidean norm, ||i?(n)i? r (n)||K, is graphed for n = (cos0,sin0) T , 0 € [0, 27r] with 
constant (p,V,p,B) using entropy scaled eigenvectors, the naturally scaled eigenvectors 
given in Sec. 1.4.3, and a slightly generalized form of the eigenvectors given in Roe and 
Balsara [RB96]. The singularities in the natural scaling are clearly seen. Once again note 
that the entropy scaled eigenvectors produce a constant result since R(n)R T (n) = u )V which 
is independent of n. 

In conclusion, using the Scaling Theorem 1.2.1, scaled eigenvectors, /? p (u), have been 
obtained for Powell’s quasi-conservative MHD equations such that 

A p (n) = Rp(n)A p (n)R 7 ; 1 (n), A 0 = I^,(n)R^(n) 
and the symmetric coefficient matrices 

A p (n)A° = Rp(n)A(n)Rj(n). 

In future sections, these results will be exploited in the construction of stabilized numerical 
discretizations. 
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Chapter 2 

Maximum Principles for Numerical 
Discretizations on Triangulated 
Domains 


One of the best known tools employed in the study of differential equations is the maximum 
principle. Any function f(x) which satisfies the inequality /" > 0 on the interval [a, 6] 
attains its maximum value at one of the endpoints of the interval. Solutions of the inequality 
/" > 0 are said to satisfy a maximum principle. Functions which satisfy a differential 
inequality in a domain ft and because of the form of the differential equation achieve a 
maximum value on the boundary dft are said to possess a maximum principle. Recall the 
maximum principle for Laplace’s equation. Let Au = u xx -\-u yy denote the Laplace operator. 
If a function u satisfies the strict inequality 

Au > 0 (2.1) 


at each point in fi, then u cannot attain its maximum at any interior point of ft. The strict 
inequality can be weakened 

Au > 0 (2.2) 


so that if a maximum value M is attained in the interior of ft then the entire function must 
be a constant with value M. Without any change in the above argument, if u satisfies the 
inequality 

Au + c\u x + c 2 u y > 0 

in then u cannot attain its maximum at an interior point. 

The second model equation of interest is the nonlinear conservation law equation: 

Ut + (/(«))*= 0, ^=a(u) (2.3) 

In the simplest setting the initial value problem is considered in which the solution is spec- 
ified along the x-axis, w(a;,0) — u 0 (. / r) in a periodic or compact supported fashion. The 
solution can be depicted in the x - t plane by a series of converging and diverging char- 
acteristic straight lines. From the solution of (2.3) Lax provides the following observation: 
the total increasing and decreasing variations of a differentiable solution between any pairs 
of characteristics are conserved [Lax73]. 


X(t T to) — /(*o), 



du(x, t) 
dx 


dx 
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Moreover in the presence of entropy satisfying discontinuities the total variation decreases 
(information is destroyed) in time. 

l(t + t Q )<I(to) (2.4) 

An equally important consequence of Lax’s observation comes from considering a mono- 
tonic solution between two non-intersecting characteristics: between pairs of characteris- 
tics, monotonic solutions remain monotonic , no new extrema are created. Also from (2.4) 
it follows that 

1. Local maxima are nonincreasing 

2. Local minima are nondecreasing 

These properties of the differential equations serve as basic design principles for numerical 
schemes which approximate them. The next section reviews the theory surrounding discrete 
matrix operators equipped with maximum principles. 

2.1 Discrete Maximum Principles for Elliptic Equations 

2.1.1 Laplace’s Equation on Structured Meshes 

Consider Laplace’s equation with Dirichlet data 

Cu = 0,x, y E fi 

u = g,x,y € dtt (2.5) 

C = A 

Prom the maximum principle property we have that 

sup|u(a;,j/)| < sup \u(x,y)\ 

:r€0 x€dQ 

For simplicity consider the unit square domain 

= {(x,y) e R 2 : 0 < x,y < 1} 

with spatial grid Xj = jAx , y^ = kAx , and J Ax = 1. Let Uj^ denote the numerical 
approximation to u(xj, i/*). It is well known that the standard second order accurate ap- 
proximation 

C&U = ^^2 Wj+hk + + U h k- i-i + Uj,k~i ~ 4 Uj 7 k] (2.6) 

exhibits a discrete maximum principle. To see this simply solve for the value at (j, k) 

U h k = ^[Uj+i,k + Uj- i 7 k + Uj,k + 1 + Uj t k- 1 ]- 
If Uj^ achieves a maximum value M in the interior then 

M + Uj-Uk + Uj,k+\ + Uj,k"i] 

which implies that 

M — Uj+ i t k = Uj—i t i z = Uj k~ i 

Repeated application of this argument for the four neighboring points yields a discrete 
maximum principle. 
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2.1.2 Monotone Matrices 


The discrete Laplacian operator -£a obtained from (2.6) is one example of a monotone 
matrix. 

Definition: A matrix M is a monotone matrix if and only if Ad'" 1 > 0 (all entries are 
nonnegative). 

Theorem 2.1.1 (Monotone Matrix) A sufficient but not necessary condition for M 
monotone is that M be an M-matrix. M-matrices have the sign pattern Mu > 0 for 
each i, Mij < 0 whenever i ^ j. In addition M must either be strictly diagonally dominant 

n 

Mu > 5Z i = 1,2 , ...,n (strict diagonal dominance) 

or else M must be irreducible and 

n 

Mu > ^2 i = 1,2, ...,n (diagonal dominance) 

3 = 

with strict inequality for at least one i. 

Proof: The proof for strictly diagonally dominant M is straightforward. Rewrite the matrix 
operator in the following form 

M = D-N, D > 0>N >0 
= [I — ND~ 1 ]D D 1 > 0 

= [I-P)D P> 0 (2.7) 

From the strict diagonal dominance of M, eigenvalues of P N I) 1 are less than unity. 
This implies that the Neumann series for [7 — P] 1 is convergent. This yields the desired 
result: 

M~ x = D~ l [I + P + P 2 + P 3 + ...} >0 (2.8) 

When M is not strictly diagonally dominant then M must be irreducible so that no per- 
mutation V exists such that 

V T MV = ^ 12 (reducibility). 

0 Mn \ 

This insures that eigenvalues of P are less than unity. Once again, the Neumann series is 
convergent and the final result follows immediately. ■ 

2.1.3 Laplace’s Equation on Unstructured Meshes 

Consider solving the Laplace equation problem (2.5) on a planar triangulation using a 
Galerkin finite element approximation with linear elemental shape functions. (Results using 
a finite volume method are identical but are not considered here.) Next pose (2.5) in vari- 
ational form. Let S h € H l be the finite-dimensional space of trial functions with bounded 
energy which satisfy the Dirichlet boundary condition on T. Similarly, let denote 

the finite-dimensional space of functions satisfying homogeneous boundary conditions. Find 
u £ S h such that for all w E V h 

[ (Vit ■ Vw) dx - 0 (2.9) 

J n 
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with 


u(x) = g(x ), x E T. 

Prom this simple equation, we have the following remarkable theorem: 


Lemma 2.1.1 The (7° linear Galerkin finite element discretization of the 2-D Laplace 
equation (2.9) is a monotone discretization if and only if the triangulation is a Delaunay 
triangulation. 

Proof: Consider a single arbitrary simplex T — simplex(xi,X 2 ,^ 3 ) and the discretization 

of (2.9) in terms of local linear shape functions N{(x) satisfying Ni(xj) = 5ij. Using these 
shape functions u(x) = Nj(x)uj, x € T and w(x) = Nj(x)wj, x € T . Inserting 
these expressions into (2.9) yields 


/, 


Vw ■ Vu dx 


3 3 

= EX>* 

i~ 1 j = 1 


Uj (VN{ • VA/j) meas(T). 


( 2 . 10 ) 


These expressions can be collected pairwise for edges surrounding a vertex. After some 
straightforward manipulation, the following global discretization formula is obtained 



where Mi denotes the set of vertices adjacent to vertex v t with weights 


( 2 . 11 ) 


Wij = (VNi • VTVj)meas(T) + (ViV' ■ VArj)meas(r') 
= — ^cotan(ajj) + cotan(aL)^ . 


(2.12) 


In this formula, a%j and a f tJ are the two angles subtending the edge e(v t ,Vj), see Fig. 2.1. 



Figure 2.1: Discretization weight geometry for the edge e(i H,vj) 

Since the discretization formula must hold for arbitrary values of Wi at interior vertices, it 
can be concluded that for all interior vertices vi 

£ (ui - Uj ) = 0. (2.13) 
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Written in this form, the discretization is monotone if all weights are nonnegative, W tJ > 0 
and Dirichlet boundary conditions are enforced. Further simplification of the edge weight 
formula is possible 


Wij 


1 

2 

1 

2 

1 

2 


^cotan(ajj) +cotan(a^)^ 


/cosjay) cosja^)\ 
y sin (ajj) sin(a^)/ 
( sin(a i j+a-j) \ 
^sin^^sin^)^ ' 


(2.14) 


Since < 7 r, a[ - < 7r, the denominator is always positive, hence nonnegativity requires 
that aij + < 7r. 



Figure 2.2; Circumcircle test for adjacent triangles, p f interior, p ,f exterior, p cocircular. 


Some trigonometry reveals that for the configuration shown in Fig. 2.2 with circumcircle 
passing through w*} the sum + ot \ • • depends on the location of p with respect to 

the circumcircle in the following way: 


a ij 

+ a 'i 3 

< 


p exterior 

a ij 

+ a ij 

> 


p interior 

a ij 

+ a ij 

= 

7T, 

p cocircular 


Also note that by considering the circumcircle passing through similar results 

would be obtained for v £. The condition of nonnegativity implies a circumcircle condition 
for all pairs of adjacent triangles whereby the circumcircle passing through either trian- 
gle cannot contain the fourth point. This is precisely the unique characterization of the 
Delaunay triangulation [Del34] which completes the proof. ■ 

Observe that from equation (2.12) that eotan(a) > 0 if a < 7t/2. Therefore a sufficient but 
not necessary condition for nonnegativity of the Laplacian weights is that all angles of the 
mesh be less than or equal to 7r/2. This is a standard result in finite element theory [CR73] 
and applies in two or more space dimensions. 

Given Lemma 2.1.1, it becomes straightforward to obtain a discrete maximum principle for 
Laplace’s equation on Delaunay triangulations using a Galerkin finite element approxima- 
tion. 
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Theorem 2.1.2 The discrete Laplacian operator obtained from the Galerkin finite element 
discretization of (2.9) with C° linear elements exhibits a discrete maximum principle for 
arbitrary point sets in two space dimensions if the triangulation of these points is a Delaunay 
triangulation. 


Proof: Prom Lemma 2.1.1, a one-to-one correspondence exists between nonnegativity of 
weights and Delaunay triangulation. Assume a Delaunay triangulation of the point set so 
that for some arbitrary interior vertex vq and all adjacent neighbors v\ we have Wo t > 0. 
Next solve for uq: 


Uq 


SieA/o WpiUj _sr^ 

^6v.wb i h 0 tl 


with 


Wi bt 


which satisfies Oi > 0 and X^eA/o G i = 1. This implies uq is a convex combination of the 
neighboring values, hence 

minuj <uq < maxuj (2.16) 


If uq attains a maximum value M then all Uj = M. Repeated application of (2.16) to 
neighboring vertices in the triangulation establishes the discrete maximum principle. ■ 


One can ask if the result concerning Delaunay triangulation and the maximum principle 
extends to three space dimensions. Unfortunately, the answer is no. This can be demon- 
strated by counterexample. The resulting formula for the three-dimensional Laplacian edge 
weight is 

^ d{v 0} Vi) 

Woi = - ]T |AR fc+i | cotan(« fc+i ). (2.17) 

0 *=i 

In this formula, A/q is the set of indices of all adjacent neighbors of vq connected by 



Figure 2.3: Set of tetrahedra sharing interior edge e(vo y vi) with local cyclic index k. 


incident edges, k a local cyclic index describing the associated vertices which form a polygon 
of degree d(vo^vi) surrounding the edge e(c 0 , u;), is the face angle between the two 

faces associated with S k+ i and which share the edge e(w*,Vfc+i) and |AR fc+ i| is the 

magnitude of the edge, see Fig. 2.3. A maximum principle is guaranteed if rill Wo t > 0. We 
now will proceed to describe a valid Delaunay triangulation with one or more Woi < 0. It 
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will suffice to consider the Delaunay triangulation of N points in which a single point Po lies 
interior to the triangulation and the remaining N - 1 points describe vertices of boundary 
faces which completely cover the convex hull of the point set. 



Figure 2.4: Subset of 3-D Delaunay Triangulation that fails to maintain nonnegativity. 


Consider a subset of the N vertices, in particular consider an interior edge incident to vq 
connecting to vi as shown in Fig. 2.4 by the dashed line segment and all neighbors adjacent 
to Vi on the hull of the point set. In this experiment the height of the interior edge, z , serves 
as a free parameter. Although it will not be proven here, the remaining N — 8 points can be 
placed without conflicting with any of the conclusions obtained for looking at the subset. 

It is known that a necessary and sufficient condition for the 3-D Delaunay triangulation 
is that the circumsphere passing through the vertices of any tetrahedron must be point free 
[Law86]; that is to say that no other point of the triangulation can lie interior to this sphere. 
Furthermore a property of locality exists so that only adjacent tetrahedra need be inspected 
for the satisfaction of the circumsphere test. For the configuration of points shown in Fig. 
2.4, convexity of the point cloud constrains z > 1 and the satisfaction of the circumsphere 
test requires that z < 2. 

1 < z < 2 (Delaunay Triangulation) 

From (2.17), Wqi > 0 if and only if z < 7/4. 

7 

1 < z < — , (Nonnegativity) 

This indicates that for 7/4 < z < 2 a valid Delaunay triangulation exists which does not 
satisfy a discrete maximum principle. In fact, the Delaunay triangulation of 400 points 
randomly distributed in the unit cube revealed that approximately 25% of the interior edge 
weights were of the wrong sign (negative). 

Keep in mind that from (2.17) the sufficient but not necessary condition for nonnega- 
tivity that all face angles be less than or equal to 7r/2. This is consistent with the known 
result from [CR73]. 

2.2 Discrete Total Variation and Maximum Principles for 
Hyperbolic Equations 

This section examines discrete total variation and maximum principles for scalar conserva- 
tion law equations. Begin by considering the nonlinear conservation law equation: 

u t + (/(«))* = u(x,0) = u 0 (x), x,t, 6 R x R + 
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which is discretized in the conservation form: 


m i+ 1 = m- ^-(h i+L - h-x) 

3 3 AX 3 ^2 3 2 7 

= H(U?_ h U?_ l+l ,...,U J+l ) (2.18) 


where /i J+ i = h(Uj-i+i, ...,Uj+i) is the numerical flux function satisfying the consistency 
condition 

h(U,U,...,U)=f(U). 

A finite-difference scheme (2.18) is said to be monotone in the sense of Harten, Hyman, and 
Lax [HHL76] if H is a monotone increasing function of each of its arguments. 

>0 V — k <i < k (HHL monotonicity) 

OU% 

This is a strong definition of monotonicity. In Crandall and Majda [CM80] it is proven that 
schemes on Cartesian grids satisfying this condition converge to the physically relevant, 
entropy satisfying solution. Kroner el ai [KRW96] has recently proven a similar result 
for monotone upwind finite volume schemes on triangulated domains. Unfortunately, HHL 
monotone schemes in conservation form are at most first order spatially accurate. Very 
few results are known concerning the convergence of high order accurate approximations. 
Johnson and Szepessy [JS90] have shown convergence to entropy solutions using streamline 
diffusion with specialized shock capturing operators. Kroner et ai [KSR95] have recently 
obtained measure-valued convergence of higher order upwind finite volume schemes for 
scalar conservation laws in several space dimensions. 

To circumvent the first order accuracy of monotone schemes, Harten introduced a weaker 
concept of monotonicity. A grid function U is called monotone if for all i 


min(Ui-i,Ui+i) < Ui < max(f/j_i, C/i+i). 

A scheme is called monotonicity preserving if monotonicity of U n + l follows from monotonic- 
ity of U n . Observe the close relationship between monotonicity preservation in time and the 
discrete maximum principle for Laplace’s equation (2.16) in space. It follows immediately 
from the definition of monotonicity preservation that 

1. Local maxima are nonincreasing 

2. Local minima are nondecreasing 

which is a property of the conservation law equation. Using this weaker form of monotonicity 
Harten [Har83a] introduced the notion of total variation diminishing schemes. Define the 
total variation in one dimension: 

oo 

TV(U) = '£ \Ui~Ui-y\. 

-oo 

A scheme is said to be total variation diminishing (TVD) if 

TV(U n+x ) < TV (U n ) 

This is a discrete analog of the total variation statement (2.4) given for the conservation 
law equation. Harten has proven that schemes which are HHL monotone are TVD and 
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schemes that are TVD are monotonicity preserving. Furthermore, it can be shown that all 
linear monotonicity preserving schemes are at most first order accurate. Thus high order 
accurate TVD schemes must necessarily be nonlinear in a fundamental way. 

To understand the basic design principles for TVD schemes, assume a one-dimensional 
periodic grid together with the following numerical scheme in abstract matrix operator form 

[/ + 0At M D]U n+l = [/—(! — 0)A t M D\U n (2.19) 


where M and M are matrices which can be nonlinear functions of the solution U. The 
matrix D denotes the difference operator 


/ 


DU = [r - E~ l ]U = 


Ui-Uj 

u 2 -Ur 

U 3 - U 2 


\ 


\Uj-Uj^J 


The scheme (2.19) represents a general family of explicit (0 = 0) and implicit (0 = 1) 
schemes with arbitrary support. More importantly, schemes written in conservative form 
can be rewritten in this form using (exact) mean value constructions. Using this notation 
an equivalent definition of the total variation in terms of the L\ norm is produced 


TV(U) = \\DUh. 


To analyze the scheme (2.19), multiply by D from the left and regroup. 

[7 + 9 At D M)D U n+1 = [/ — (! — 0)A t D M]D U n 


or in symbolic form 


CD U n+1 =nDU n D U n+l = C l HD U n 

where invertibility of C has been assumed. This invertibility will be guaranteed from the 
diagonal dominance required below. Next take the L\ norm of both sides and apply matrix- 
vector inequalities. 


TV(U n+1 ) = ||DU n+1 ||i < ||/r 1 ft||i||DT / n || 1 

= j|£- - 1 ll\\iTV{U n ) (2.20) 

This reveals the sufficient condition that \\C 1 1?- 1 J i < 1. Recall that the L\ norm of a 
matrix is obtained by summing the absolute value of elements of columns of the matrix 
and choosing the column whose sum is largest. Furthermore, we have the usual matrix 
inequality ||£" l 7£||i < ||£ _l ||i||72.||i so that sufficient TVD conditions are that ||£ l ||i < 1 
and \\1Z\\ i < 1. These simple inequalities are enough to recover the TVD criteria of previous 
investigators, see [Har83b, JL86]. 

Theorem 2.2.1 A sufficient condition for the scheme (2.19) to be TVD with 0 = 0 is that 
TZ be bounded with 1Z > 0 (all elements are nonnegative). 
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Proof: Consider the explicit operator TZ = I — At D M and multiply it from the left by 
the summation vector s 1 = [1, 1, ..., 1], It is clear that s T D = 0 so that s T TZ — s 1 (columns 
sum to unity). Because the Li norm of TZ is the maximum of the sum of absolute values of 
elements in columns of 1Z, the stated theorem is proven. ■ 

Next consider the implicit scheme with 0=1 and sufficient conditions for \\C “MU < i- 
Prom the previous development, one way to do this would be to show that C is monotone, 
i.e. CT X > 0 with columns that sum to unity. 

Theorem 2.2.2 A sufficient condition for the scheme (2.19) to be TVD with 0 = 1 is that 
C be an M-type monotone matrix, i.e. diagonally dominant with positive diagonal entries 
and negative off-diagonal entries. 
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To obtain Harten’s TVD criteria for the implicit scheme, we need only require that this 
operator be an M-matrix to obtain the following conditions as did Harten 


D Un £ 0 


D- + 1/ 2 > o. 


2,2.1 Maximum Principles and Monotonicity Preserving Schemes on Mul- 
tidimensional Structured Meshes 

Unfortunately, two motivations suggest a further weakening the concept of monotonicity. 
The first motivation concerns a negative result by Goodman and Le Veque [GV85] that 
conservative TVD schemes on Cartesian meshes in two space dimensions are first order 
accurate. The second motivation is the apparent difficulty in extending the TVD concept to 
arbitrary unstructured meshes. The first motivation inspired Spekreijse [Spe87] to consider 
a new class of monotonicity preserving schemes based on positivity of coefficients. Consider 
the following conservation law equation in two space dimensions 


u t + (/(«))* + (9(u))y = 0 . 

Next construct a discretization of (2.21) on a logically rectangular mesh 






At 


= A U^> k - u ” k) + a j hM >■* - 


( 2 . 21 ) 


( 2 . 22 ) 


with nonlinear coefficients 



B ± 

3 


1 

2 


,k 


1,*, 

lifc , 


u? >k ,u?+ hk , 


...) 

...) 


Theorem 2.2.3 The scheme (2.22) exhibits a discrete maximum principle at steady-state 
if all coefficients are uniformly bounded and nonnegative 




> 0 B 


J ± 2’ k 


> 0. 


Furthermore, the scheme (2.22) is monotonicity preserving under a CFL-like condition if 

A . 1 

At < min — r x r- 


Proof: The first task is to prove a discrete maximum principle at steady-state by solving 
for the value at (j y k). 


Uj* 


i.* + Bj,k±±Uj,k± 0 
52±( A j±±,h + 

± 1 


(2.23) 


25 



with the constraints 


a i-{M + a J+h k + + ^> k +b k ~ 1 

and a j± i k > 0, and (i- k± ± > 0. Prom positivity of coefficients and convexity of (2.23) it 
follows that 

f^fedbi) < Uj y k 5: max(f/j : ti j / c , Uj }k ±\)* (2.24) 

If U jtk attains a maximum value M at (j y k) then 


M — Uj~ — Uj+i t k — Ujfi-i — U jMl . 


Repeated application of (2.24) to neighboring mesh points establishes the maximum prin- 
ciple. 

Next it is straightforward to establish a CFL-like condition for monotonicity preservation 
in time by again seeking positivity of coefficients and a convex local mapping from U n to 
f/ n+1 . 


Kt' 


(l - A V?, k + A« T.(A 1±i , k UJ ±lx + 

+ Y^( a j±i,k^J± l,k + fij,k±jU” k±1 ) 


(2.25) 


with the derivable constraints 


7j,k + a j-±,k + a j+$,k + + Pjjk+i* ~ 1 

and aj ± i k > 0, and /?. k±k > 0. To show that (2.25) is a local convex mapping, it suffices 
to satisfy the CFL-like condition for nonnegativity of 7 ^: 


At < min — - —c 

-VjkZ±(AZ 




(2.26) 


If (2.26) is satisfied then monotonicity preservation in time follows immediately: 
min (U? ±l , k ,U? MU U? tk ) < < max(U? ±lik ,U” k±l ,U? >k ). 


2.2.2 Maximum Principles and Monotonicity Preserving Schemes on Un- 
structured Meshes 

This section examines the maximum principle theory for conservation laws on unstructured 
meshes. Specifically, our primary attention focuses on Godunov-like upwind finite volume 
schemes [God 59 ] utilizing solution reconstruction and evolution. Some early maximum prin- 
ciple results for upwind finite volume schemes can be found in [DD 88 ] [BJ89] [Bar91]. Note 
that many of these results were subsequently used in implementations of the discontinuous 
Galerkin method as well, see for example Bey [Bey91]. Also note that the present analysis 
differs from maximum principle theory based on the “upwind triangle” scheme developed 
by Desideri and Dervieux [DD 88 ], Arminjon and Dervieux [AD93]. 
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Consider the integral conservation law form of (2.21) for some domain f2 comprised of 
nonoverlapping control volumes, flj, such that SI — Lift* and ft* n f ^ = 0 , i ^ j. Next, 
impose the integral conservation law statement on each control volume 

^ [ udQ+ [ ( Fn)dT = 0 (2.27) 

ot Jn i Jan, 

where F(u ) = / (u)i + g(u)j. The situation is depicted for a control volume fto in Fig. 2.5. 
For two- and three-dimensional triangulations, several control volume choices are available: 



Figure 2.5: Local control volume configuration for unstructured mesh. 


the triangles themselves, Voronoi duals, median duals, etc. Although the actual choice of 
control volume tessellation is very important, the monotonicity analysis contained in the 
remainder of this section is largely independent of this choice. Consequently, a generic 
control volume fto with neighboring control volumes ft*, i E A/o, as shown in Fig. 2.5, is 
sufficient for the present analysis. In the following example, the solution data is assumed 
constant in each control volume. This simplifies the analysis considerably. The second 
example addresses the more general situation utilizing high order data reconstruction. 

Example: Analysis of an Upwind Finite Volume Scheme with Piecewise Con- 
stant Reconstruction 

In this example, assume that the solution data u(x , y)i in each control volume ft* is constant 
with value equal to the integral average value, i.e. 

u(x^y)i = u* = -j- / u dft, Vf It E ft 

Ai J n t 

where A* is the area of ft*. Next, define the unit exterior normal vector n o* for the control 
volume boundary separating fto and f 2*. It is also useful to define a normal vector no * 
which is scaled by the length (area in 3-D) of that portion of the control volume boundary 
separating fto and ft*. Finally, to simplify the exposition, define 

/(«;n) - F(u) ■ n 

and assume the existence of a mean value linearization such that 

f(v;n) -f(u;n) = df (u,v\n)(v - u). (2.28) 
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(2.29) 


Using this notation, construct the following upwind scheme 

^t(Aouo) = ~ 

<tt i<zN 0 


with 

/i(w 0) Wi;n 0 t) = - (f(u 0 ;noi) + f(u ; n oi )) - -|d/(w 0 ,^;n 0 «)|(wj - w 0 ) 

In Barth and Jespersen [BJ89], we proved a maximum principle and monotonicity preser- 
vation of the scheme (2.29) for scalar advection. 


Theorem 2.2.4 The upwind algorithm (2.29) with piecewise constant solution data ex- 
hibits a discrete maximum principle for arbitrary unstructured meshes and is monotonicity 
preserving under the CFL-like condition: 


At < min 




vajeQ EteMj 4f' («*>«* 5 “*) 


Proof: Consider the control volume surrounding no as shown in Fig. 2.5. Recall that the 
flux function was constructed using a mean value linearization such that 

/(«»; no;) - f(u 0 ; noi) = df(u 0 ,uf,n 0 i) ( Tii -uo) 

This permits regrouping terms into the following form: 

^-(u 0 A 0 ) = Y 1 /(«o;noi) - 5Z df (uo,Ui;n 0 i)(ui -u 0 ) 

1 teM o ieAfo 


where (*) = (•)~ l ~ + (•) and |(*)| = («) + — (•) . For any closed control volume, it follows that 

Y2 /(«o;n 0 j) = 0. 
ieAfo 


Combining the remaining terms yields a final form for analysis 

^-(uqAq) = - y' d/-(uo,Uj;noi)(ui - u 0 ). 
i^Afo 


(2.30) 


To verify a maximum principle at steady-state, set the time term to zero and solve for uq. 


Uq 


(^o^ijnQj) Uj __ 

EieVo d f ~ («0 > ; n 0t ) a,Ul 


with Y^ie A/o a * = 1 and a* > 0. Since uq is a convex combination of all neighbors 


min Ui < uq < ma x u*. (2.31) 

ieMo ieAf 0 v ’ 


If t«o takes on a maximum value M in the interior, then n, = M,V i G A/o- Repeated 
application of (2.31) to neighboring control volumes establishes the maximum principle. 
For Euler explicit time stepping, 


£ 

dt 


(?t 0 A 0 ) « v4 0 


u? +1 


0 


0 


At 
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a CFL-like condition is obtained for monotonicity preservation. Inserting this expression 
into (2.30) yields 


At 


a 


o +1 = «o - ~r H d f («0,<;n 0 »)« - «o) 

ieA/'o 

= ao«o + Y aiU i ■ 

ieJVo 


(2.32) 


It should be clear that coefficients in (2.32) sum to unity. To prove monotonicity preserva- 
tion, it is sufficient to show nonnegativity of coefficients. Clearly, c*j > 0 V i > 0, hence, 
monotonicity preservation is achieved if 

«o = i + ^ Y > o. 

710 i€A/o 

This implies monotonicity preservation in time under the CFL-like condition 


At < min 


-A* 


VOjen df not) 


Example: Analysis of High Order Accurate Upwind Advection Schemes Using 
Arbitrary Order Reconstruction [Bar94] 

In this example, high order accurate upwind schemes on unstructured meshes are consid- 
ered. The technique used here is to show a maximum principle for the cell averages. The 
solution algorithm is a relatively standard procedure for extensions of Godunov’s scheme in 
Eulerian coordinates, see for example [God59, vL79, CW84, HOEC87]. The basic idea in 
Godunov’s method is to treat the integral control volume averages, u, as the basic unknowns. 
Using information from the control volume averages, k — th order piecewise polynomials are 
reconstructed in each control volume ftf. 

U k (x,y)i= Y, a {m,n) P (m,n ) (* ~ *c, 1/ ~ Vc ) 

7re-f n<k 

where P( mj7J ) ( x — x Cy y — y c ) — (x — x c ) m (y — y c ) n and (x c , y c ) is the control volume centroid. 
The process of reconstruction amounts to finding the polynomial coefficients, Near 

steep gradients and discontinuities, these polynomial coefficients maybe altered based on 
monotonicity arguments. Because the reconstructed polynomials vary discontinuously from 
control volume to control volume, a unique value of the solution does not exist at control 
volume interfaces. This non-uniqueness is resolved via exact or approximate solutions of the 
Riemann problem. In practice, this is accomplished by supplanting the true flux function 
with a numerical flux function which produces a single unique flux given two solution states. 
Once the flux integral is carried out (either exactly or by numerical quadrature), the control 
volume average of the solution can be evolved in time. In most cases, standard techniques 
for integrating ODE equations are used for the time evolution, i.e. Euler implicit, Euler 
explicit, Runge-Kutta. The result of the evolution process is a new collection of control 
volume averages. The process can then be repeated. The process can be summarized in the 
following steps: 
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(1) Reconstruction in Each Control Volume: Given integral solution averages in all 
Clj, reconstruct a k —th order piecewise polynomial U k (x,y)i in each fi* for use in equation 
(2.27). In faithful implementations of Godunov’s method, cell averages of the solution data 


[ U k (x,y) t dQ = (uA) t 
Jvti 


are preserved during the reconstruction process. For solutions containing discontinuities 
and/or steep gradients, monotonicity enforcement may be required. 

(2) Flux Evaluation on Each Edge: Supplant the true flux by a numerical flux function. 
Given two solution states the numerical flux function returns a single unique flux. Using 
the notation of the previous section, define f(u ; n) — ( F(u ) • n) so that 



dr ^ f h(u L ,u R -,n)dr. 

JdQi 


Consider each control volume boundary to be a collection of polygonal edges (or dual 
edges) from the mesh. Along each edge (or dual edge), perform a high order accurate 
flux quadrature. When the reconstruction polynomial is piecewise linear, single (midpoint) 
quadrature is usually employed on both structured and unstructured meshes 

h(U L ,U R ;n)dV « E h(U L ,U R -,n) ij 



where U L and U R are solution values evaluated at the midpoint of control volume edges as 
shown in Fig. 2.5. When multi-point quadrature formulas are employed, they are assumed 
to be of the form: 

f(s)ds = fitq) 

q£Q 

with w q > 0 and £ [0, 1]. Let the multi-point quadrature formulas be represented by the 
augmented notation 



n)dr« E E 

jeAfiqeQ 


(3) Evolution in Each Control Volume: Collect flux contributions in each control 
volume and evolve in time using any time stepping scheme, i.e. Euler explicit, Euler implicit, 
Runge-Kutta, etc. The result of this step is once again control volume averages and the 
process can be repeated. 

In the present analysis, the reconstruction polynomials U k (x, y)i in each $1* are given. 
The result of the analysis will be conditions or constraints on the reconstruction so that 
a maximum principle involving cell averages can be obtained. The topic of reconstruction 
and implementation of the constraints determined by this analysis will be examined in 
a later section. Using this notation, the following upwind scheme is constructed for the 
configuration in Fig. 2.5. 


4(^o«o) = - E E Wqh(U L ,U R ;n)ou, 

ieM, q€Q 


(2.33) 
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with a numerical flux function obtained from (2.28). 


h{U‘\ U R ; noa) = \ (j f(U L ; n) + f(U R ; n)) Q . 

- ^\df(U L ,U R -,n)\oi(U R -U L ) 0i (2.34) 

To analyze the scheme, recall that the flux function was constructed using a mean value 
linearization such that 

f(U R ; n) - /(t/ L ; n) = df(U l % U R ; u){U R - U L ). 


This permits regrouping terms into the following form: 


j t {u,A 0 ) = - E J2 w i(f( uL -’ A )+ d r( uL ’U R -,n)(U R -U L )) o . . (2.35) 

*6A/o<?eQ lq 

Rewrite the first term in the sum using a mean value construction 

E ^ w qf(^o',n)oi q + w q (df(u 0 ,U L -,n)(U L -uo)) n . . 

ie-VoqeQ V /0tq 

The first term vanishes when summed over a closed volume so that (2.35) reduces to 


Ao) = - E E (df(uo,U L -,n)(U L -uo) + dr{U L ,U R -,n){U R -U L )) 

ieAfoqeQ J 

By introducing difference ratios, the scheme can be written in the following form: 


0 iq 


= “EE™? ( d f («o,tf L ;n)¥) ( u i~ u o) 

ieJVoqeQ iq 

- E E w 9(^ + ^ ,[/t; "^)oi ^ o_5t ) 

ieAfoqeQ xq 

- EE w q (< <V~(U L , U R ; n)0) o . (Ui - u 0 ) 


ieAfo q€Q 


(2.36) 


with 


,t, U olq ~ «o ^ - n 0 ^ v oi - 1 j 0i 

$ 0*9 — — =r— , $ 0 a 9 = — r— , © 0*9 = — r r: — ■ 

It ■ 1 1 ^ ■-> I ~ <11, ^ .» J . r, i _ 


u R n - uL 

$0*9 = ©0*9 = 

Ui — Uo Uo — Uk Uj — u 0 

In this equation, the k subscript refers to some as yet unspecified index value, k € Mo- 


Theorem 2.2.5 The generalized Godunov scheme with arbitrary order reconstruction (2.33) 
exhibits a discrete maximum principle at steady-state if the following three conditions are 
fulfilled: 

$U 9 > °> ®jiq > 0 ®jiq > 0 Y?, 4 , i € Mj (2.37) 

as defined by (2.36). Furthermore, the scheme is monotonicity preserving under a CFL-like 
condition if 


At < min 


-At 


v ii 


en EieAO Z qeQ - d f + * + 4re) .. 9 
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Proof: Assume that (2.37) holds. Define df oiq = w q df(uo, U L \ n)oiq and similarly dfo iq = 
w q dj{U Ij > U R ] n)ozg- Setting the time term to zero and solving for no yields 


_ g i^o^Q - (<r +*■ 

E q€Q 9- <T e)^ 

= ^ " atiUi. (2.38) 

ieMi o 


Examining the individual coefficients, it is clear that a i = 1 and e*i > 0,Vi. Thus a 

convex local mapping exists and 


minimi < z^o < maxt/*. (2.39) 

ieAo i6A/*o 


If no takes on a maximum value M in the interior, then n* — M,V i € A/o- Repeated 
application of (2.39) to neighboring control volumes establishes the maximum principle. 

To establish monotonicity preservation in time, consider Euler explicit time-stepping 
scheme. 


s (uo/io) 


^4o 


TK +l 




At 


Inserting this formula into (2.36) yields 


tt«+i - 


u n 


s-(r E E^*®.! 5 ?-^) 

' l0 teMi q eQ 

7 E E ^ 1 0 ! 9 ^ 0 i ?(«0 — u fc ) 

tEE d foi q ®Otq(u 7 i - BS) 

^ i6A/o 

+ E E 


(2.40) 


with c*o + Xwt.v'n = 1 and a, > 0, i > 0. A locally convex mapping in time from U n to 
U n+l is achieved when ao > 0. This assures monotonicity in time. Some algebra reveals 
the following formula for «o 


oc o = l + 


At 

Ao 


EE Or 

ieAfa qeQ 


<H-df + <f> + 



0 iq 


From this the CFL-like condition for monotonicity preservation is obtained 

a j ^ . ~Ao 

At < min 7= — — 

v ih€ " E,ee (4f «p - 4f * + <T e)^ 

so that 

?nin(«”,iio) < < majc(uf ,Ug) 

leNo ie No 

Applying this result, to all control volumes establishes monotonicity preservation. ■ 

Without specifying the actual type of reconstruction, we have the following simple lemma 
concerning the ratios appearing in (2.36): 
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Lemma 2.2.1 Assume 'I'oig > 0 as defined in (2.36). A sufficient condition for <&ot <7 > 0 is 
that the reconstructed polynomial reduce to the cell average value, U k (x,y ) o = uo, at local 
cell average extrema, i.e. whenever 


max u i < iin < minti,. 
jeM 0 J J 


Proof: Consider an arbitrary control volume flj adjacent to fio- Assume that rtj > uq. 
The stated assumption, 'Joiq > 0 implies that U^ q > u 0 . Consequently, $o J9 < 0 if and 
only if uq is less than all adjacent neighbors, hence a local minimum. Following a similar 
argument, uq is a local maximum when Ui < 7I () and 4>oiq < 0. ■ 
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Chapter 3 

Upwind Finite Volume Schemes 


This chapter examines upwind finite volumes schemes for scalar and system conservation 
laws. The basic tasks in the upwind finite volume approach have already been presented: 
reconstruction, flux evaluation, and evolution. By far, the most difficult task in this process 
is the reconstruction step. 

3.1 Reconstruction Schemes for Upwind Finite Volume Schemes 

In the following paragraphs, the design criteria for general reconstruction operators with 
fixed stencil is reviewed. The reader is referred to the papers by Abgrall [Abg94, Abg95], 
Vankeirsbilck [Van93] and Michell [Mic94] for a discussion of ENO and adaptive stencil 
reconstruction schemes. 

The reconstruction operator serves as a finite-dimensional (possibly pseudo) inverse of 
the cell-averaging operator A whose j-th component A j computes the cell average of the 
solution in ftj. 



In addition, the following properties are usually imposed on the reconstruction: 

(1) Conservation of the mean: Simply stated, given cell averages n, we require that all 
polynomial reconstructions u k have the correct cell average. 

if u k = R k u then u = A u k 

This means that R fc is a right inverse of the averaging operator A. 

AR fc = I 

Conservation of the mean has an important implication. Unlike finite-element schemes, 
Godunov schemes have a diagonal mass matrix. 

(2) k-exactness: A reconstruction operator R fc is A '-exact if R*A reconstructs polynomials 
of degree k or less exactly. 

if u E Vk and u — A u, then u k — R k u = u 

In other words, R fc is a left-inverse of A restricted to the space of polynomials of degree at 
most k . 

R fc A = I 

Vk 
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This insures that exact solutions contained in Vk are in fact solutions of the discrete equa- 
tions. For sufficiently smooth solutions, the property of ^-exactness also insures that when 
piecewise polynomials are evaluated at control volume boundaries, the difference between 
solution states diminishes with increasing A; at a rate proportional to /i fc+l were h is a 
maximum diameter of the two control volumes. Figure 3.1 shows a global quartic polyno- 
mial u e V 4 which has been averaged in each interval. This same figure shows linear and 



Figure 3.1: Cell averaging of quartic polynomial (left), linear reconstruction (center) and 
quadratic reconstruction (right). 


quadratic reconstruction given interval averages. The small jumps in the piecewise polyno- 
mials at interval boundaries would decrease even more for cubics and vanish altogether for 
quartic reconstruction. Property (1) requires that the area under each piecewise polynomial 
is exactly equal to the cell average. 

One of the most important observations concerning linear reconstruction is that one can 
dispense with the notion of cell averages as unknowns by reinterpreting the unknowns as 
pointwise values of the solution sampled at the centroid ( midpoint in 1-D) of the control 
volume . This well known result greatly simplifies schemes based on linear reconstruction. 
The lineal' reconstruction in each interval shown in Fig. 3.1 was obtained by a simple 
central-difference formula given pointwise values of the solution at the midpoint of each 
interval. Note that for steady-state computations, conservation of the mean in the data 
reconstruction is not necessary. The implication of violating this conservation is that a 
nondiagonal mass matrix appears in the time integral. Since time derivatives vanish at 
steady-state, the effect of this mass matrix vanishes at steady-state. The reconstruction 
schemes presented below assume that solution variables are placed at the vertices of the 
mesh, which may not be at the precise centroid. 

3.1.1 Green-Gauss Reconstruction 

Let Do denote the set of all triangles incident to some vertex vq and the exact integral 
relation 

[ Vudtt= [ undT . (3.1) 

JDo JdDo 
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It is not difficult to show [BJ89] that given function values at vertices of a triangulation, a 
discretization of this formula can be constructed which is exact whenever u varies linearly: 


(V«)„ fl = -j- 53 + tio)noi. 

A ° ieM o 1 


(3.2) 


In this formula not = J'f d n for any path which connects triangle centroids adjacent to the 
edge e(vo, Vj) and Ao is the area of the nonoverlapping dual regions formed by this choice of 
path integration. Two typical choices are the median and centroid duals as shown below. 
This approximation extends naturally to three dimensions. The formula (3.2) suggests a 

4 

3 

Mesh 

Median Dual 
Centroid Dual 


Figure 3.2: Local mesh with centroid and median duals. 



natural computer implementation using the edge data structure. Assume that the normals 
n ij for all edges e(vi,Vj) have been precomputed with the convention that the normal vector 
points from v\ to Vj . An edge implementation of (3.2) can be performed in the following 
way: 

For k = 1 ,n(e) / Loop through edges of mesh 
ji = e~ 1 (k i 1) / Pointer to edge origin 
j 2 = e _1 (A:,2) ! Pointer to edge destination 
uav = -h u(j 2))/2 ! Gather 

ux(jy) -f = nonnx{k) • uav ! Scatter 
uxfa) — = normx(k) ■ uav 
u v{h) T = normy(k) ■ uav 
u y{h) = normy(k) ■ uav 
End for 

For j — l,n(i?) / Loop through vertices 

ux(j) = ux(j) /area(j) ! Scale by area 
uy(j) = uy(j)/area(j) 

End for 

It can be shown that the use of edge formulas for the computation of vertex gradients is 
asymptotically optimal in terms of work done. 

3.1.2 Linear Least-Squares Reconstruction 

To derive this reconstruction technique, consider a vertex vq and suppose that the solution 
varies linearly over the support, of adjacent neighbors of the mesh. In this case, the change 
in vertex values of the solution along an edge e(v^vo) can be calculated by 


(Vtx)o • (r, - r 0 ) = Ui -u 0 . 
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This equation represents the scaled projection of the gradient along the edge e(v;,uo)- A 
similar equation could be written for all incident edges subject to an arbitrary weighting 
factor. The result is the following matrix equation, shown here in two dimensions: 

'uqAaq W\Ay\ 

_w n Ax n w n Ay n 

or in symbolic form C Vu = f where 

C = [L l L 2 ) 



in two dimensions. Exact calculation of gradients for linearly varying u is guaranteed if any 
two row vectors u>i(rj - ro) span all of 2 space. This implies linear independence of L\ and 
L 2 . The system can then be solved via a Gram-Schmidt process, i.e., 


A 

A. 


[L i 


L 2 ) = 


1 0 
0 1 


The row vectors V are given by 


Vi 


V 2 


I22L1 — I12L2 
I11I22 ~ l\ 2 
I11L2 — hjLi 
^11^22 — ^12 


( 5 . 3 . 4 ) 


( 3 . 3 ) 


with lij = (Li ■ Lj). 

Note that reconstruction of N independent variables in R d implies +dN inner 

product sums. Since only dN of these sums involves the solution variables themselves, 
the remaining sums could be precalculated and stored in computer memory. This makes 
the present scheme competitive with the Green-Gauss reconstruction. Using the edge data 
structure, the calculation of inner product sums can be calculated for arbitrary combinations 
of polyhedral cells. In all cases linear functions are reconstructed exactly. This technique 
is best illustrated by example: 

For k = l,n(e) / Loop through edges of mesh 
ji = e~ 1 (k, 1 ) / Pointer to edge origin 
j2 = e~ l (k, 2 ) / Pointer to edge destination 
dx = w(k) ■ (x(j2) — x(j 1)) / Weighted Ar 
dy = w(k) ■ (y(h) - 2 /(ji» Weighted Ay 
ln(jt) = hi(ji) + dx ■ dx ! I u orig sum 

ln(j2) = hi(h) + dx-dx ! l u dest sum 

h2(ji) = h 2(31) + dx ■ dy ! l l2 orig sum 

h 2 (h) = U2U2) + dx • dy ! l l2 dest sum 

du = w(k) ■ (u(./2) — u(ji)) / Weighted An 
Ifi (ji ) + = dx ■ du !Lif sum 
>f\(h) + = dx ■ du 
lh(j\) + = dy - du !L 2 f sum 
If 2(h) + = dy- du 
Endfor 

For j = 1 , n(v) ! Loop through vertices 
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det = l n {j) I22U) ~ lit 

ux(j) = (I22U) ■ ~ h2 ■ If 2) I det 

uy(j) = (ln(j)-lf2(j)-h2lfi)/det 

End for 

This formulation provides freedom in the choice of weighting coefficients, wi. These 
weighting coefficients can be a function of the geometry and/or solution. Classical approx- 
imations in one dimension can be recovered by choosing geometrical weights of the form 
Wi = l./| Arj — Aro|* for values of t — 0, 1 , 2. 

3.1.3 Monotonicity Enforcement 

When solution discontinuities and steep gradients as present, additional steps must be 
taken to prevent oscillations from developing in the numerical solution. One way to do this 
was pioneered by van Leer [vL79] in the late 1970’s. His basic idea was to enforce strict 
monotonicity in the reconstruction. Monotonicity in this context means that the value of the 
reconstructed polynomial does not exceed the minimum and maximum of neighboring cell 
averages. The final reconstruction must guarantee that no new extrema have been created. 
When a new extremum is produced, the slope of the reconstruction in that interval is 
reduced until monotonicity is restored. This implies that at a local minimum or maximum 
in the cell averaged data the slope in 1-D is always reduced to zero, see for example Fig. 
3.3. Theorem 2.2.5 provides sufficient conditions for a discrete maximum principle in the 



Figure 3.3: Linear data reconstruction with monotone limiting. 


cell averages using arbitrary order reconstruction on general control volumes. Consider the 
control volume interface separating flo and as shown in Fig. 2.5. From Theorem 2.2.5, 
a maximum principle is guaranteed if for all quadrature points on the interface separating 
Qo and $ 1, 

'Foi > 0, $o; >0 0 O i > 0. 

Lemma 2.2.1 states that <f»o £ > 0 is always satisfied if the monotonicity enforcement algo- 
rithm reduces to piecewise constant at local extremum, i.e. when 

max m, < Tin < min u , . 
jeM 0 3 - ~ jtM'o 3 
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Assume that this property holds, monotonicity reduces to the following two conditions at 
all quadrature points: 

°< fej? « 

0< tHt M < 3 - 4 ) 

The second inequality appearing in (3.4) requires that the difference in the extrapolated 
states at a cell interface must be of the same sign as the difference in the cell average 
values. For example in Fig. 3.4(a) this condition is violated but can be remedied either by 
a symmetric reduction of slopes or by replacing the larger slope by the minimum value of 
the two slopes. Observe that in one space dimension the net effect of the slope limiting in 



Figure 3.4: (a) Reconstruction profile with increased variation violating monotonicity con- 
straints. (b) Profile after modification to satisfy monotonicity constraints. 

the 7'econstruction process is to ensure that the total variation of the reconstructed function 
does not exceed the total vanation of the cell averaged data. 

In Barth and Jespersen [BJ89], a simple recipe was proposed for slope limiting of lin- 
early reconstructed data on arbitrary unstructured meshes. Consider writing the linearly 
reconstructed data in the following form for 

U{x,y) o = uo + Vu 0 • (r - r 0 ). 

Now consider a “limited” form of this piecewise linear distribution. 

u (*, y ) o = U 0 + $0 V«o • (r - r 0 ) 


u’ 0 nm = min (u 0 ,Ui) 
t€N o 

Uo nax = max(«o,«») 


and require that 


«r‘ < U(x,y ) o < tC 


when evaluated at the quadrature points used in the flux integral computation. For each 
quadrature point location in the flux integral, compute the extrapolated state Ufc and 
determine the smallest 4>o so that 




inin(l, o L ~ ), if V < - it 0 > 0 

u 0t “0 

..min tt 

<min(l ^ ^ < 0 

u 0t u ° 

k 1 if U& - u 0 = 0 
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This limiter function automatically satisfies Lemma 2.2.1. Condition (a) from (3.4) is local 
to the control volume and can be enforced by a further reduction in slope. In practice this 
step is sometimes omitted. Condition (b) is enforced using the procedure shown in Fig. 3.4. 

Extensive numerical testing has shown that this limiter can noticeably degrade the over- 
all accuracy of computations, especially for flows on coarse meshes. In addition the limiter 
behaves poorly when the solution is nearly constant unless additional heuristic parameters 
are added. This has prompted other researchers (c.f. Venkatakrishnan [Ven93]) to proposed 
alternative '‘smooth’ 1 limiter functions, but no serious attempt is made to appeal to the rig- 
ors of maximum principle theory. The design of accurate limiters satisfying the maximum 
principle constraints is still an open problem in this area. 

When the above procedures are combined with the flux function given earlier (2.34), 

h(U\U a ;n) = \(f(U L -,n)+f(U R ;n)) 

- l -\df(U L ,U R ^)\(u R -U L ) (3.5) 

the resulting scheme has good shock resolving characteristics. To demonstrate this, we 
consider the scalar nonlinear hyperbolic problem suggested by Struijs, Deconinck, et al 
[Str89]. The equation is a multidimensional form of Burger’s equation. 

■u t + (u 2 / 2) x + u y = 0 

This equation is solved in a square region [0,1.5] x [0,1.5] with boundary conditions: 
u(x,0) — 1.5 — 2x , x < 1, u(x,0) = —.5, x > 1, u(0,y) — 1.5, and u(1.5,y) — —.5. 
Figure 3.5 shows carpet plots and contours of the solution on regular and irregular meshes. 
The exact solution to this problem consists of converging straightline characteristics which 



Figure 3.5: Carpet plot of Burger’s equation solution on (a) regular mesh, (b) irregular 
mesh. 


eventually form a shock which propagates to the upper boundary. The carpet plots indicate 
that the numerical solution on both meshes is monotone. Even so, most people would prefer 
the solution on the regular mesh. This is an unavoidable consequence of irregular meshes. 
The only remedy appears to be mesh adaptation. 

Next, consider a test problem which solves the two-dimensional scalar advection equation 


u t + {yu) x - (xu) y = 0 
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or equivalently 


ut + A • Vu — 0, A — (y, —x) T 

on a grid centered about the origin, see Fig. 3.6. Discontinuous inflow data is specified 
along an interior cut line, u(x,0) = 1 for —.6 < x < —.3 and u(x,0) = 0, otherwise. The 
exact solution is a solid body rotation of the cut line data throughout the domain. The 




Figure 3.6: Grid for the circular advection problem (left) and solution contours obtained 
using piecewise constant reconstruction (right). 



Figure 3.7: Solution contours, piecewise linear reconstruction (left) and piecewise quadratic 
reconstruction (right). 

discontinuities admitted by this equation are similar to the linear contact and slip line so- 
lutions admitted by the Euler equations. Figure 3.6 also shows solution contours obtained 
using piecewise constant reconstruction. The discontinuities are quickly smeared by the 
computation. Figure 3.7 displays solution contours obtained using piecewise linear and a 
piecewise quadratic reconstruction technique discussed in [Bar93, Bar94]. The improvement 
from piecewise constant reconstruction to piecewise linear is quite dramatic. The improve- 
ment from piecewise linear to piecewise quadratic also looks impressive. The width of the 
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discontinuities is substantially reduced with little observable grid dependence. Note how- 
ever, that the quadratic approximation used for this computation has roughly quadrupled 
the number of solution unknowns because of the use of 6-noded triangles. 


3.2 Numerical Solution of the Euler and Navier-Stokes Equa- 
tions Using Upwind Finite Volume Approximation 

The section discusses the extension and application of the scalar finite volume scheme to 
the Euler and Navier-Stokes equations. One advantage of the finite volume method is the 
relative ease in which this can be done. 


3.2.1 Extension of Scalar Advection Schemes to Systems of Equations 

The extension of the scalar advection schemes to the Euler (and Navier-Stokes) equations 
requires three rather minor modifications: 


1. Vector Flux Function . The scalar flux function is replaced by a vector flux function. 
In the present work, the mean value linearization due to Roe [Roe81] is used. The 
form of this vector flux function is identical to the scalar flux function (2.34), i.e. 

h(u ft , u L ; n) = i (f (u R , n) + f(u L ; n)) 

- ||A(u fl ,u L ;n)| (u ft -u L ) (3.6) 

where f(u; n) = F(u) • n, and A = df/du is the flux Jacobian. 

2. Componentwise limiting. The solution variables are reconstructed componentwise. 
In principle, any set of variables can be used in the reconstruction (primitive vari- 
ables, entropy variables, etc.). Note that conservation of the mean can make certain 
variable combinations more difficult to implement than others because of the non- 
linearities that may be introduced. The simplest choice is obviously the conserved 
variables themselves. When conservation of the mean is not important (steady-state 
calculations), we typically use primitive variables in the reconstruction step. 


3. 


Weak Boundary Conditions . Boundary conditions for inviscid flow at solid surfaces 
are enforced weakly. For solid wall boundary edges, the flux is calculated with V * n 
set identically to zero 


f(u;n) 


( 0 \ 

n x p 


. n yP . 

Vo/ 


Boundary conditions at far field boundaries are also done weakly. Define the charac- 
teristic projectors of the flux Jacobian A in the following way: 


P* = ±[I±siffi(A)]. 


At far field boundary edges the fluxes are assumed to be of the form: 


f(u”; n) = (F(u" roj ) - n) 
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where u£ r0J = P+ u n + P Uoo and u represents a vector of prescribed far field 
solution values. At first glance, prescribing the entire vector is an overspecification 
of boundary conditions. Fortunately the characteristic projectors remove or ignore 
certain combinations of data so that the correct number of conditions are specified at 
inflow and outflow. 

3.2-2 Example: Supersonic Oblique Shock Reflections 

In this example, two supersonic streams (M=2.50 and M=2.31) sire introduced at the left 
boundary. These streams interact producing a pattern of supersonic shock reflections down 
the length of the converging channel. The grid is a subdivided 15x52 mesh with perturbed 
coordinates as shown in Fig. 3.8. This same figure also shows Mach contours for the nurner- 



Figure 3.8: (a) Channel grid, 780 vertices (left), (b) Mach contours, piecewise constant 
reconstruction (right). 


ical solution obtained using piecewise constant reconstruction. As expected, the piecewise 
constant reconstruction scheme severely smears the shock system while the scheme based 
on a linear solution reconstruction shown in Fig. 3.9 performs very well. The piecewise 
quadratic approximation, shows some improvement in shock wave thickness although the 
improvement is not dramatic given the increased number of unknowns used in this partic- 
ular scheme [Bar93]. The number of unknowns required for the quadratic approximation is 
roughly four times the number required for the piecewise linear scheme. This lack of dra- 
matic improvement is not a surprising result since the solution has large regions of constant 
flow which do not benefit greatly from the quadratic approximation. At solution disconti- 
nuities the quadratic scheme reduces to a low order approximation which again negates the 
benefit of the quadratic reconstruction. To provide a fair comparison, Fig. 3.10 shows the 



Figure 3.9: Mach contours, (a) Piecewise linear reconstruction (left), (b) Piecewise 

quadratic reconstruction (right). 


same mesh adaptively refined. The number of mesh points has roughly doubled. A numeri- 
cal solution was then obtained using linear reconstruction. The results are very comparable 
to the calculation performed using quadratic reconstruction while requiring less that 10% 
as much computational time. 
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Figure 3.10: (a) Adapted channel grid, 1675 vertices, (left), (b) Mach contours, piecewise 
linear reconstruction (right). 


3.2,3 Example: Transonic Airfoil Flow 

Figure 3.11 shows a simple Steiner triangulation and the resulting solution obtained with a 
linear reconstruction scheme for transonic Euler flow = .80, a = 1.25°) over a NACA 
0012 airfoil section. 



Even though the grid is very coarse with only 3155 vertices, the upper surface shock 
is captured cleanly with a profile that extends over two cells of the mesh. Clearly, the 
power of the unstructured grid method is the ability to locally adapt the mesh to resolve 
flow features. Figure 3.12 shows an adaptively refined mesh and solution for the same flow. 
The flow features in Fig. 3.12 are clearly defined with a weak lower surface shock now 
visible. Figure 3.13 shows the surface pressure coefficient distribution on the airfoil. The 
discontinuities are monotonically captured by the scheme. 

3.2.4 Example: Navier-Stokes Flow with Turbulence 

As a last finite volume example, compressible Navier-Stokes flow is computed about a 
multiple-component airfoil geometry. In addition to the basic Navier-Stokes equations, 
the effects of turbulence on the mean flow equations are modeled using an eddy viscosity 
turbulence transport model. In a report with Baldwin [BB90], we proposed a single equation 
turbulence transport model with the specific application to unstructured meshes in mind. 
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Figure 3.12: (a) Solution adaptive triangulation of airfoil, 6917 vertices, (b) Mach number 
solution contours on adapted airfoil. 



Figure 3.13: Comparison of C p distributions on initial and adapted meshes. 


This model was subsequently modified by Spalart and Allmaras [SA92] to improve the 
predictive capability of the model for wakes and shear-layers as well as to simplify the 
model’s dependence on distance to solid walls. In the present computations, the Spalart 
model is solved in a form fully coupled to the Navier-Stokes equations. The one-equation 
model for the viscosity-like parameter u is written 

= 1 [ v • (( u + £)V£) +C 62 (VP) 2 ] - Cwxjn +C b iSv. (3.7) 

In the Spalart model the kinematic eddy viscosity is given by is t = uf v \ and requires the 
following closure functions and constants 

— i 

■X , , ISIS Y 

S = aj| + ~2~nfv 2) /ul = ■ 3 , 3 

X'* + Si 
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fv 2 

<7 


l+Xfv i’ 
r + c u , 2 (r 6 - r) 


5« 2 <f 2 


with w the fluid vorticity, d the distance to the closest surface, and the constants c&i = 
0.1355, c 62 = 0.622, q,! = 7.1, = 3.24, 2 = 0.3, c„,3 = 2.0, k = .41, <7 = 2./3.. 

The model also includes an optional term for simulating transition to turbulence. Using 
this model, viscous flow with turbulence is computed about the multiple-element airfoil 
geometry. This geometry has been triangulated using the Steiner triangulation algorithm 



Figure 3.14: Multi-element airfoil mesh (left) and solution isomach contours (right), M 0 0 = 
0.2, a = 16.0°, Re = 9.0 million. 



Figure 3.15: Comparison of computation and experiment (left), solution convergence history 
(right) 

described in [Bar95a], see Figure 3.14. The relatively coarse mesh contains approximately 
22,000 vertices with cells near the airfoil surface attaining aspect ratios greater than 1000:1. 
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This example provides a demanding test case for CFD algorithms. The experimental flow 
conditions are M = .20, a = 16°, and a Reynolds number of 9 million. Experimental 
results are given in [VDMG92]. Even though the wake passing over the main element is 
not well resolved, the surface pressure coefficient shown in Figure 3.15 agrees quite well 
with experiment. The spatially discretized flow equations are solved using an exact Newton 
interation as described in [Bar95b] . The eventual quadratic convergence of Newton’s method 
is shown in Fig. 3.15. 
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Chapter 4 

Simplified GLS in Symmetric Form 


This chapter returns to the symmetrization topic of Chapter 1 and presents a new simplified 
variant of the Galerkin least-squares (GLS) scheme in symmetric form. By combining the 
Scaling Theorem 1.2.1 of Chapter 1 with a simple congruence approximation given in the 
next section, a substancial simplification of this scheme is possible. 

4.1 Congruence Approximation 

Lemma 4.1.1 (Congruence Approximation) Let A £ R nx " be an arbitrary diagonal- 
izable matrix and B £ R Ttxn a right symmetrizer of A. If R £ R nx ” denotes the right 
eigenvector matrix scaled under Theorem 1.2.1 such that 

A — RAR~ l , and B = RR T , 

then the following inertially equivalent approximate decomposition exists 

(dB) ± « RA ± R t 


satisfying 


AB = (AB) + + (AB)~ = RA + R T + RA~ R T . 


Proof: Let U £ R nxn denote the unitary matrix diagonalizing AB 

ab = unu T = RAR t , UU T = /, RR t = B. 

Left and right multiplication by U T and U respectively yields 

a = (u t r)a(u t r) t . 


Prom Sylvestor’s law of inertia 


Inertia(Q) = Inertia(A) 


so consequently 

Inertia(i i^) — Inertia^*) = Inertia (^{U 1 R)A ± (U 1 R) 1 ^ . 
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Left and right multiplication by U and U T respectively 

UtfU 7 ' = (AB) ± “ RA ± R t . 

Finally, adding positive and negative components reveals 

R( A+ + A ~)R T = RAR t = A B. 


4.2 Simplified Galerkin Least-Squares in Symmetric Form 

Let 7" =]t n ,t n+1 [ denote the nth time interval and il the spatial domain composed of 
nonoverlapping elements T u G = UTj, 1} Pi Tj = 0, i / j. Next define the trial and weighting 
finite element spaces 

S h = {v' l |v /l £ (c°(Q x /"))"*, v*| T x/» £ (v k (T X /"))”*, q(v) = g D on V D x 7 n } 

v h = {w fc |w fc £ (c°(n x /"))"*, w A | T x/«» e (r k (T x 7 n )) m ,q , (w) = 0 on T x 7 n } 

where v denotes the entropy variables for the system. The Galerkin least-squares method 
is defined by the following bilinear forms [HFM86]: 

Find E V h such that for all E S h 


B(v\ v/ h ) ga i + B(v h ,w h )i s + B(v h , w h ) bc = 0 


with 


B(v,w) gal = f j ( u(v) ■ W( f*(v)w ja;i dCl, dt 
Ji n Jet 

+ [ (w(r +1 )u(v(P +1 )) - w(t”)u(v(<” ))) dn 

J n 

f ^2 f (A°w t + A‘A°w iXl j r (vl 0 v t + A l A°v ;Xi ^ dYldt 
Jln re q Jt 

[ [ w h(v,g F ;n) dr dt 

Ji n J r F 


B(v,w)i s 

B(v,w)bc 


where 


h(v..,v+,n) = i(f(u(v._);n) + f(u(v+); n)) - ^|A(u(v); n)|(u(v+) - u(v_)). 


4.2.1 A Simplified Least-Squares Operator in Symmetric Form 


In the implementation of the Galerkin least-squares method, the most difficult computa- 
tional aspect of the scheme is the calculation of the least-squares r matrix. In the paper by 
Hughes and Mallet [HM86], they proposed the following general form for r on a mapped, 
x i-» £, anisotropic element 


r 


p 


= \B\ 


-l 

p » 


\B\ P = 



i/p 

, B l = (d^/dxjjAJA 0 


(4.1) 
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where A 3 and A 0 are the matrices given in (1.5). Historically, the value p — 2 has been 
used in implementations of the Galerkin least-squares method. After manipulation, t 2 can 
be written in the form 



This necessitates the calculation of a matrix square root. Hughes advocates the use of 
the Cayley-Hamilton theorem for this computation. Unfortunately, the Cayley-Hamilton 
technique becomes unwieldy for matrices of dimension larger than 4 or 5. In light of the 
Scaling Theorem 1.2.1, it is useful to revisit the derivation of r with p = 1. Define n* = 
V&/IV6I and A(n) = ntA 1 so that (4.1) can be rewritten as 

T! = l^lr 1 = (|V6IM(n i M°|)' 1 . 

Observe that for the Euler and MHD equations considered in Chapter 1, the eigenvectors 
of A(n)A° are not easily obtained. Hence, the computation of |A(ni)A°| seems nontrivial 
as well. Lemma 4.1.1 suggests the inertially equivalent approximation 

|A(n).4 0 | « R(n)\A{n)\R T (n) 

using the entropy scaled eigenvectors R(n) of /l(n). From this, the following approximation 
for r follows 

ti « (|V&| ft(nj)|A(nj|# T (nj)) . 

This represents a substancial simplification of the r matrix calculation. 

4.2.2 Example: MHD Flow for Perturbed Prandtl-Meyer Fan 

In this example, the MHD equations are solved in the square domain ft € [—1/2, 1/2] 2 using 
the simplified GLS scheme. The inflow data consists of a Mach 1.55 Prandtl-Meyer fan with 
origin located at (a; = —.881 , y = .47147). A velocity aligned magnetic field, B = .05V, is 
introduced at the boundary. Figure 4.1 shows the coarsest mesh used for the computation 
and Mach number contours computed from the simplified GLS scheme. Figure 4.2 shows 
contours of the x-compoiient of the magnetic field using linear and quadratic elements on 
the coarsest mesh. Solutions were then obtained on a sequence of meshes with decreasing 
mesh spacing and the L 2 norm of V « B measured. Figure 4.3 graphs the convergence of 
V*B for linear and quadratic approximations. The graphs show optimal convergence rates, 
0(h ) for linear and 0(h 2 ) for quadratic elements. 
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Figure 4.1: Coarsest mesh (400 vertices) and Mach number contours for perturbed Prandtl- 
Meyer fan. 




Figure 4.2: B x component of magnetic field computed on coarsest mesh. Galerkin least- 
squares with linear elements (left) and quadratic elements (right). 



Figure 4.3: Spatial Convergence of V • B for the Galerkin least-squares scheme using linear 
and quadratic elements. 
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